function [E,freq,velocity] = phaseShiftTransform(data,Fs,distance,RangeFreq,ParaVelocity,nu,dataMode)
%   计算面波的频散能量图，具备相位加权叠加功能
%   输入参数:
%       data：地震数据，大小为[nt*nx]，其中nt是时间采样点数，nx是接收点数量。
%       Fs: 信号采样率。
%       distance: 接收点到震源的距离数组，大小为[nx]。
%       RangeFreq: 色散能量中的频率范围（Hz），例如 [5 100]。
%       ParaVelocity: 相速度（m/s）的参数，第一个是最小值，第二个是最大值，第三个是间隔，
%                     例如 [50 800 1]。
%       nu：默认取0，使用则通常取值为1。
%       dataMode：归一化方式。0表示列归一化，1则是全局归一化。
%   输出参数:
%       E: 归一化的频散能量矩阵。
%       freq: 频散能量图中的频率向量。
%       velocity: 相速度向量。

% @author:夜剑听雨
% E-mail:2530595378@qq.com
% 2025.08.29

% 设置默认参数
if nargin < 6 
    nu = 0; % 默认不启用相位加权
    dataMode = 0; % 默认列归一化
end

if nargin < 7 
    dataMode = 0; % 默认列归一化
end
% 计算FFT参数
[nt,~] = size(data);  % 获取地震记录时间长度

% 设定频率分辨率控制填充长度
delta_f = 0.25;   % 设置频率分辨率为0.25Hz
N_fft = Fs/delta_f;  % 采样率是偶数，可确保填充长度也是偶数

if nt > N_fft
    data = data(1:N_fft, :);
elseif nt < N_fft
    data = [data; zeros(N_fft - nt, size(data, 2))];
end

% 沿时间维度进行快速傅里叶变换
fft_data = fft(data, N_fft, 1);

% 提取正频率部分（包含零频率）
positive_freq_data = fft_data(1:N_fft/2+1, :);

% 构建正频率索引
freq_ind = (0:N_fft/2) * delta_f;  % 正频率的频率值向量，从0开始

% 确定RangeFreq(1)-RangeFreq(2)对应的索引范围
index_min = find(freq_ind >= RangeFreq(1), 1, 'first');  % 找到大于等于最小频率的第一个索引
index_max = find(freq_ind <= RangeFreq(2), 1, 'last');  % 找到小于等于最大频率的最后一个索引

% 截取RangeFreq(1)-RangeFreq(2)的目标数据
obj_data = positive_freq_data(index_min:index_max, :);
% 对数据频谱作归一化处理
obj_data = obj_data ./ abs(obj_data);

% 对应的频率值
freq = freq_ind(index_min:index_max);

% 根据给定参数构建相速度向量
velocity = ParaVelocity(2):-ParaVelocity(3):ParaVelocity(1);

% 初始化频散能量矩阵
lf = length(freq);       % 频率的数量
lv = length(velocity);   % 相速度的数量
E = zeros(lv, lf);


% 遍历频率和速度，计算相位加权能量
for j = 1:lf
    for i = 1:lv
        k0 = 2 * pi * freq(j) / velocity(i);
        phase_shift = exp(1i * k0 * distance);
        phase_shifted_signals = phase_shift .* obj_data(j,:);

        % 计算归一化复数信号（幅度为1）
        normalized_signals = phase_shifted_signals ./ abs(phase_shifted_signals);
        normalized_signals(isnan(normalized_signals)) = 0; % 处理除零异常

        % 计算全局相位一致性因子（公式5）
        coherence = abs(mean(normalized_signals));  % 各道归一化信号的复数平均模长
        W = coherence ^ nu;  % 应用幂指数

        % 复数叠加并应用全局权重（公式6）
        S = sum(phase_shifted_signals);  % 线性叠加
        E(i,j) = abs(S) * W;  % 相位加权后的能量


    end

end

% 归一化频散能量图
if dataMode==0
    for j=1:lf
        % 归一化处理
        E(:,j) = abs(E(:,j)) ./ max(abs(E(:,j)));
    end
else
    E = abs(E) ./ max(abs(E));
end
% 处理能量矩阵中的NaN或Inf值
E(isnan(E)) = 0.0; % 将NaN值替换为0
E(isinf(E)) = 0.0; % 将Inf值替换为0
end

%%   参考文献:
%       Park CB, Miller RD, Xia J. Imaging dispersion curves of surface waves
%       on multi-channel record[C]. Technical Program with Biographies SEG,
%       68th Annual Meeting, New Orleans, L A., 1998, 1377-1380.
%       Cheng, F., Xia, J., Zhang, K., Zhou, C., & Ajo-Franklin, J. B. (2021).
%       Phase-weighted slant stacking for surface wave dispersion measurement.
%       Geophysical Journal International, 226(1), 256–269. https://doi.org/10.1093/gji/ggab101

